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ABSTRACT 

The most massive neutron stars constrain the behavior of ultra-dense matter, with larger masses possible only 
for increasingly stiff equations of state. Here, we present evidence that the black widow pulsar, PSR B 1957+20, 
has a high mass. We took spectra of its strongly irradiated companion and found an observed radial-velocity 
amplitude of K a \, s = 324 ± 3kms _1 . Correcting this for the fact that, due to the irradiation, the center of light 
lies inward relative to the center of mass, we infer a true radial-velocity amplitude of K2 = 353 ±4kms _I and 
a mass ratio q = Mpsr/M2 = 69.2 ± 0.8. Combined with the inclination i = 65 ± 2 deg inferred from models of 
the lightcurve, our best-fit pulsar mass is Mpsr = 2.40 ± 0. 12M©. We discuss possible systematic uncertainties, 
in particular in the lightcurve modeling. Taking an upper limit of i < 85 deg based on the absence of radio 
eclipses at high frequency, combined with a conservative lower-limit to the motion of the center of mass, 
K2 > 343 kms" 1 (q > 67.3), we infer a lower limit to the pulsar mass of Mpsr > 1 .66M Q . 
Subject headings: pulsars: individual (PSR B 1957+20) — stars: neutron 



1, INTRODUCTION 

One of the outstanding problems in physics is the behav- 
ior of matter at extreme densities. This behavior is modelled 
using quantum-chromodynamics calculations, but these can- 
not yet determine reliably the densities at which, e.g., meson 
condensation and the hadron to quark-gluon phase transition 
occur. At densities slightly above nuclear and high tempera- 
ture, models can be tested with heavy-nuclei collision experi- 
ments. For higher densities and low temperatures, only com- 
pari son with neutron-star param eters is possible (for a review, 
e.g.. lLattimer & Prakashll2007l) . 

The different models lead to different equations of state, 
which predict different mass-radius relations for neutron stars. 
Unfortunately, most attempts at observational tests have been 
frustrated by susceptibility to systematic errors and modeling 
uncertainties. The most robust tests have involved measure- 
ments of extrema. For i nstance, the fastest measured spin 
period, 1.4 ms (Ter 5ad; iHessels et al.1 120061) . excludes the 
stiffest equations of state, for which neutron stars would be 
too large to spin so fast. 

A problem with measurements of the extrema, is that 
whether they occur in Nature depends not only on whether 
they are allowed physically, but also whether they are ex- 
pected astron o mical ly. From models of stellar evolution, 
Timmes et al. (1996) find neutron-star mass distributions at 
birth with two narrow peaks, at 1.3 and 1.8M Q , containing 
remnants of stars with initial masses smaller and larger than ~ 
19M©, respectively. In binaries, however, where much of the 
stellar envelope is removed during the evolution, they expect 
only to form neutron stars in the lower mass bin. The above 
may explain why until recently most accurate masses were all 
close to 1.4M©: most were measured for binaries containing 
pulsars with neutron star companions, where the preceding 
evolution predicts rel atively little mass has been accreted (for 
a review, Stairs 2004). The one exception was the X-ray bi- 
nary Vela X-l, for which a higher mass of 1.86±0.16M© 



was inferred dBarziv et al.ll200lHOuaintrell et al.ll2003l) . Such 
a high mass would imply, e.g., that meson condensation is not 
important. However, the large uncertainty prevented a defini- 
tive conclusion. 

However, neutron stars can become more massive after 
birth by accretion. Accretion also leads to "recycling" of radio 
pulsars: it increases their spin frequency and decreases their 
magnetic fields. Most recycled pulsars are accompanied by 
low-mass, <0.2M© companions. For these systems to form 
in the age of the Universe, the companion must originally 
have been a star of >0.8M Q . Thus, the companion lost al- 
most all of its mass, and some of it should have landed on the 
neutron star, leading to a concomitant increase in the neutron- 
star mass. Initial tests, however, were suggestive but not con- 
clusive: for pulsars with low-mass white dwarf companions, 
masses around 1.4— 1.7 M© were foun d (Ivan Kerkwiik et al 
1996; [Jacobv et al.1 120051 iBassa et alj 120061: IVerbiest et al 



2008), somewhat larger than typical, but not yet very con- 
straining. Presumably, in these systems much of the mass 
actually left the binary. 

More recently, higher masses were found in different types 
of binaries, st arting with a number of pulsar binaries in glob- 
ular clusters ( Ransometal] |2005l l^eireeral] 1200833' . In 
these cases, however, the masses rely on observations of pe- 
riastron advance, which is assumed to be due to general rela- 
tivistic effects only (rather than classical ones such as due to 
rotationally and tidally induced quadrupoles), and statistical 
arguments that the inclinations are unlikely to be very low. 
Thus, it was still possible to doubt that very massive neutron 
stars could exist. 

While we were writing and revising this work, however, 
such doubts disappeared, with accurate mass determinations 
for PSR J1614-2230 (1 .97 ± 0.04M Q , iDemorest et aTJl20iOh 
and PSR J 1903+0327 (1 .67 ±O.O2M ), both relying on mea- 
surements of Shapiro delay, which is not easily mimicked by 
other processes. These masses exclude many of the soft equa- 
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tions of state, such as Kaon condensation as envisaged by 
iBrown & Beth3 ( fT994h . 

Intriguingly, both of the above systems do not have 
low-mass white-dwarf companions like most binary pul- 
sars, but rather a more massive, carbon-oxygen white dwarf 
(PSR J 16 14-2230) and a solar-mass main sequence star 
(PSR J1903+0327). Thus, both systems also had different 
evolutionary histo ries (puzzling for PSR J1903+0327; see 
iFreire et al. 2010). In this sense, our approach is similar, in 
that we specifically target another group of binary pulsars with 
different properties and evolutionary histories, and, therefore, 
perhaps different masses, the so-called "black widow pul- 
sars." 

In black-widow systems, a millisecond pulsar is accom- 
panied by a low-mass, few 0.01 Mq companion, which 
is bloated and strongly irradiated by the pulsar, lead- 
ing to outflows strong enough to eclipse the pulsar sig- 
nal for significant fractions of the orbit. The irradiation 
causes strong heating on the side of the companion fac- 
ing the pulsar, and, as a result, stro ng orbital brightness 
variations of the optica l counterparts dKulkarnietalJll988t 
Ivan Paradij set alJ 119881: IStappers et al] 119961) . From de- 
tailed modellin g of the lightcurves, the inclinations can 
be constrained (Callana netaT] 119951: IStappers et all 119991: 
iRevnolds et al.l 120071) . which, when combined with velocity 
information, can be used to derive masses. 

Here, we present a radial-velocity study of PSR B 1957+20 
dFruchter et al.l 119881) . the proto-type black-widow system, 
which has the brightest and best-studied counterpart. In §12 
we describe our observations and data reduction, and in §[3]we 
constrain the properties of the companion from the spectra, 
determine radial velocities, and fit an observed radial-velocity 
amplitude. In §|4] we discuss the available constraints on the 
radius and inclination, and in §0 the corrections we need to 
make because our velocity amplitude is that of the center of 
light, which, because of the irradiation, is shifted towards the 
pulsar relative to the center of mass. We present our final con- 
straints on the masses in §|6] and discuss how these may be 
made more secure in the future. 

2. OBSERVATIONS & REDUCTION 

We obtained pilot observations at the Keck telescope of the 
companion of PSR B1957+20 on the night of 15 June 2007, 
and a larger set of spectra on the nights of 4 and 5 August 
2008 (see TableHJ. For all observations, the seeing was good, 
ranging from 0.6 to ~ 1", but the sky was not photometric. 
For relative flux calibration, we obtaine d spectra of a num- 
ber of spectrophotometric standards from Bo hlin et all ( 1995) 
(see below). We obtained exposures of internal flat fields and 
Hg/Kr/Ar and Ne/Ar arc spectra interspersed with the obser- 
vations. 

The spectra were obtained using the two-armed Low 
Resolution Imaging Spectrometer (LRIS, lOke et al] 119951: 
iMcCarthv et al.lll998l) . We employed the atmospheric disper- 
sion corrector, used a 0"7 slit, set to position angle 35° to 
cover both the companion and a nearby star (at ~ l."3, here- 
after the contaminato£|), and split the light with a dichroic 
at 6800 A. In the blue arm, we used the 600 line mm" 1 grism, 
blazed at 4000 A, which covers 3 100-5600 A at a resolution 
AA~3.2Aor Av~220kms _1 (forthe0"7 slit). The detector 
is a mosaic of two Marconi CCDs, each with 4096 x 2048 pix- 

1 The pulsar's proper motion has increased the separation over that quoted 
in earlier publications. 
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Log of Observations and Velocity Measurements 
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2007 Jun 15 


11:41 


930 


54266.48677 


0.55 


0.65 


0.74 


419± 


6 




11:58 


930 


54266.49844 


0.58 


0.81 


0.90 


405 ± 


5 




14:40 


930 


54266.61091 


0.87 


1.13 


1.08 


-97 ± 


6 




14:54 


640 


54266.62052 


0.90 


0.94 


1.07 


-139 ±22 


2008 Aug 4 06:26 


1500 


54682.26797 


0.08 


0.18 


0.29 


-179 ±13 




10:12 


1800 


54682.42481 


0.49 


0.41 


0.50 


429 ± 


6 




10:42 


1800 


54682.44608 


0.54 


0.64 


0.69 


422 ± 


5 




11:16 


1800 


54682.46951 


0.60 


0.90 


0.90 


385 ± 


5 




11:47 


1800 


54682.49124 


0.66 


1.09 


1.05 


298 ± 


5 




12:12 


1000 


54682.50830 


0.71 


1.23 


1.18 


212± 


6 




13:06 


900 


54682.54576 


0.80 


1.31 


1.18 


22 ± 


7 




13:24 


1200 


54682.55845 


0.84 


1.18 


1.11 


-48 ± 


5 




13:50 


1800 


54682.57640 


0.88 


1.00 


1.02 


-125 ± 


6 


2008 Aug 5 


08:37 


1800 


54683.35900 


0.93 


0.75 


0.87 


-184 ± 


8 




09:08 


1800 


54683.38036 


0.99 


0.57 


0.63 


-213 ± 


9 



NOTE. — Col. (1): Date of the observation. Col. (2): Start time. Col. (3): 
Integration time. Col. (4): Mid-exposure, barycentric Modified Julian Date. 
Col. (5): Phase using epoch of as cending node T = MJD 48196.0635242 and 
orbital period P = 33001.91484 I Arzoumanian, Fruchte r7&" Tavlorl 1 994 we 
ignored orbital-period derivatives, see text). Col. (6): Flux ratio relative to 
the contaminator in the 5 100-5300 A range. Col. (7): Flux ratio in the 8450- 
8650 A range. Col. (8): Radial velocity relative to the contaminator inferred 
from the blue spectra (we estimate a barycentric velocity of the contaminator 
of-25.6±1.3kms-'). 

els of 15 fj,m on the side (0."135 on the sky), which we binned 
by two in the dispersion direction (the only direction it can be 
binned). On the red side, we used the 1200 line mm" 1 grating, 
blazed at 7500 A, set to cover 7600-8900 A at AA = 2.1 A or 
Av ~ 75 km s" 1 . Here, the detector was a Tektronix CCD with 
2048 x 2048 pixels of 24/im on the side (0"215 on the sky), 
which we read out unbinned. 

For the reduction, we used the Munich Image Data Analysis 
System ESO-MIDAS, and routines running in the MIDAS en- 
vironment. For all images, we subtracted bias as determined 
from the overscan regions. For the blue images, we subse- 
quently corrected for small-scale variations in efficiency by 
dividing by a spatially averaged flat field, normalized using a 
third-degree polynomial, and with the bluest, poorly exposed 
part shortward of 4000 A replaced by unity. For the red im- 
ages, we simply divided by the flat field, normalized using a 
bi-linear fit. 

The spectra and their uncertainties were extracted by fit- 
ting, at each dispersion position, the sum of three stellar pro- 
files (for the companion, contaminator, and another star at 
~3."2 on the other side of the target) and a constant sky. 
For the profiles, we used Moffat functions of the form P = 
A/(l+(x/w) 2 ) 5 , with power 5 = 5 for the blue images and 6 = 6 
for the red ones. At each dispersion position, only the three 
amplitudes A and the sky were fitted; the other parameters 
were determined globally (with the central position and width 
w allowed to vary quadratically with dispersion position, and 
the relative positions fitted as constants). The fits to the im- 
ages were generally good for the pulsar fields, with reduced 
X 2 near unity, and somewhat poorer for the higher signal- 
to-noise images of the spectrophotometric standards (where 
slight mismatches between the Moffat function and the true 
point spread function are more apparent). 

Wavelength calibration was done using arc spectra. For the 
blue arm, we used well-exposed images of the arc lamps taken 
at the start of the night to define an overall solution, which re- 
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quired a fourth-order polynomial to give an adequate disper- 
sion solution, with root-mean-square residuals of 0.14 A for 
24 lines. Next, we found offsets relative to this solution from 
less well-exposed arc frames taken interspersed between the 
observations. For the red arm, the individual arc frames were 
well exposed, and a third-degree polynomial sufficed to give 
solutions with residuals of ^0.04 A (for typically 31 lines). 

For flux calibration, we first corrected all spectra approxi- 
mately for atmospheric e xtinction using a curve made by com- 
bining the CFHT values dBeland. Boulade. & Davidgell 19881) 
shortward of 5200 A with the better sampled La Silla values 
longward of 5200 A (ESO users manual 1993; see also iTud 
119771) . Next, for the blue spectra, we calculated response 
curves by comparing our observed s pectra for Feige 110 with 
the calibrated STIS spectra given bv lBohlin et aTM995l) : we 
slightly smoothed our spectra to match the STIS resolution, 
divided the two, and smoothly interpolated the ratio. For the 
2007 red spectra, we proceeded similarly, except that we com- 
pared our observation of BD +28 42 1 1 with a third-degre e 
polynomial fit to the STIS spectrum of lBohlin et all d 19951) . 
since the STIS spectrum was relatively noisy and the intrin- 
sic spectrum should be smooth. Furthermore, we took care 
to fit the ratio spectrum avoiding telluric absorption, so that 
we could use the deviations - scaled with airmass - to correct 
to first order for telluric absorption in our target spectra. For 
the 2008 red spectra, our procedure was similar, except that 
we compared our observed spectrum of Hz 43A with the cal- 
ibrated model spectrum (which is again smo oth over the rel- 
evant wavelength range; Bohlin et al.| [l995). Since the con- 
ditions were not photometric, the above gives good relative 
fluxes over the observed wavelength range, but only approxi- 
mate absolute fluxes. To place the companion spectra on the 
same flux scale, we scaled all sets of spectra to the same aver- 
age contaminator flux in the 5100-5300 A and 8450-8650 A 
ranges (for blue and red, respectively). After this scaling, in- 
tegrating over B filter band passes, we find that our inferred 
companion magnitudes are roughl y consistent with the ones 
observed bv lRevnolds et all (120071) . 

3. SPECTRA AND VELOCITIES 

In Fig.Q] we show averaged blue and red spectra for phases 
within 0.15 from superior conjunction of the companion, 
when the irradiated side is seen face on, and for phases be- 
tween 0. 15 and 0.3 from superior conjunction, when the view 
is more sideways. 

Below, we first determine the spectral type for the two 
phases. While this will be only an "average" spectral type 
over the surface, it gives a way to determine the influence of 
irradiation. We will find that the spectrum appears rather nor- 
mal, suggesting that the irradiated energy is dissipated well 
below the photosphere (in contrast to the case of, e.g., the 
pre-cataclysmic variable NN Ser, where the companion is ir- 
radiated with ultra violet light from a hot white dwarf; see 
iParsons et al]|2010l) . 

Given the normal appearance of the spectra, we use them to 
help determine reddening and distance, and to constrain the 
surface gravity and radius of the companion. In §|5] we will 
find that these provide additional evidence that the companion 
is close to filling its Roche lobe. Next, we describe how we 
determined the velocities from the individual spectra, and use 
these to fit an orbit. We discuss the required corrections to 
this orbit in §|5] 



3.1. Spectral Type, Reddening, and Distance 

We compared our blue spectra with classification spectra 
shown in the on-line atlas by R. O. Gray0 We find that the 
contaminator has spectral type Gl-2 V, while that of the com- 
panion is slightly earlier at its brightest phase (F9-G0, as seen 
from, e.g., the weaker G band near 4300 A), and similar when 
seen from the side (about Gl). Its luminosity class is slightly 
higher, intermediate between IV and III, as can be seen from 
the stronger Sr II A4077 line, and also supported by the lower 
strength of the Fe I A4046 line relative to the Mn I A4030 
line. Consistently, in Fig. Q] one sees that (3 Hyi, which has 
spectral type G2IV, is intermediate in luminosity class be- 
tween the contaminator and the companion. For f3 Hyi, the 
mean density and radius have been measured using astero- 
seismology and interf erometry, leading to a surface gravity 
logg = 3.952 ± 0.005 (iNorth et al.ll2007l) . For the companion 
of PSR B 1957+20, we thus infer logg < 4. 

The blue spectra of the companion show emission cores 
in the Ca II H and K lines, suggestive of an active chromo- 
sphere. Furthermore, H/3 seems rather weak, possibly due to 
being filled in by poorly subtrac ted emission from the Balm er- 
dominated bow shock nebula (Ku lkarni & Hesterlll988t the 
contaminator might also be affected by this). Finally, the 
Mg lb triplet is a bit weaker than expected. 

In the red spectra, we find that the spectrum of the contam- 
inator is as expected for a Gl-2 V star, but that the spectra 
of the companion cannot be classified as easily. In particular, 
the Ca II IR triplet nearly absent, perhaps being filled in by 
chromospheric emission. Furthermore, the Na I AA8183,8195 
lines are much stronger than seen in the contaminator: com- 
bined equivalenth widths of ^0.9 and 1.5 A for the face-on 
and more sideways view of the companion, respectively, and 
~0.7 A for the contaminator. Probably, this reflects that in the 
red, a large contribution to the light arises from cooler regions 
of the companion, which will have much stronger Na I ab- 
sorption (e.g.. IZhoull 1 99 ll finds equivalent widths of ^0.7 A 
and ~ 1.5 A for early G and early M dwarfs, respectively). 
The strength of the Na I line also shows that the surface grav- 
ity of the companion is similar to that of a main-sequence or 
sub-giant star; for giants, the Na I equivalent width does not 
exceed 1 A (IZhoull 19911) for any temperature. 

We can use the spectral information to constrain the red- 
dening to the source. Starting with the contaminator, from 
archival Hubble Space Telescope o bservations wit h WFPC2, 
which we analyzed with HSTphot (Dolphin 2000), we mea- 
sure m F675w = 19.968 ± 0.01 and m F814w = 19.458 ± 0.015. 
Using the transformations of Holtz man et alj (119951) . this cor- 
responds to R = 20.06 ± 0.03 and /= 19.40 ± 0.02. Since a 
G1-2 V star has (R-I) = 0.34 ±0.02 and M R ~ 4.3 dCoxl 
2000), we infer an extinction Ay = 1.42 ± 0.1 8 mag and dis- 
tance of ~ 8 kpc (where we use the extinction curve of 
Schl egel. Finkbeiner. & D"avis|[l998l) . 

The reddening to the companion to PSR B 1957+20 should 
be similar to that to the contaminator, since its observed colors 
at maximum are about as much bluer (from our blue spectra, 
A(B-V) = -0.07 ± 0.02) as the intrinsic color difference ex- 
pected from the spectral type (A(B-V) = -0.05 ± 0.02). The 
reddening can be compared with the run of reddening with 
di stance measured using re d-clump stars. Using the technique 
of iDurant & van Kerkwijkl (120061) . we find that the reddening 

2 |http: //nedwww . ipac . caltech . edu/leve!5/Gray/Gray_contents .ht 
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FIG. 1 . — Spectra of the companion to PSR B 1957+20. We show averages both around superior conjunction (blue trace), when the illuminated side is most in 
view, and nearer quadratu re (red trace). Also shown are the spectrum of the nearby "contaminator" (green trace, offset by 5 units), as well as a UVES spectrum 
of the G2IV star f} Hyi (Bagnulo et al. 2003) (black trace; convolved to the same resolution, scaled to match the brightness of the conjunction spectrum, and 
offset by 10 units; the red part of that spectrum is not shown, as it has residual ripples and does not cover the Ca II triplet). The right-hand enlargement shows the 
Ca II H and K lines, which have cores in emission for the companion, and lines of Sr II at 4077 A, Fe I at 4046 A, and Mn I at 4030 A, which we use to estimate 
the luminosity class and surface gravity. 



is Ay — 0.5 at d ~ 1.5kpc, increases first slowly and then 
more sharply to Ay — 1.2 at d ~ 2kpc, and remains constant 
thereafter. Given that the contaminator is certainly well be- 
yond 2 kpc, and that the reddening of PSR B 1957+20 and the 
contaminator are similar, we conclude that PSR B 1957+20 is 
at a distance d > 2 kpc. This is consistent with the distance 
in the range of 1.5 to 2.5 kpc inferred from the pulsar dis- 
persion measure of 29cm~ 3 pc using models of the Galactic 
electr on distribution (Tay lor & Cord es 1993; Cordes & Laziol 
120021) . From the observations of ICallanan et al.1 (11995b . the 
companion at maximum brightness is about 0.4 mag brighter 
than the contaminator in the R band, which, using the R-band 
magnitude aboveEl implies R = 19.7. This corresponds to an 
absolute magnitude Mr ~ 7.2-5 log d% and an effective radius 
R ~ 0.25d 2 R Q (where d 2 = of /2 kpc). 

3.2. Radial Velocities 

We determined velocities by fitting the flux-calibrated blue 
spectra of both the pulsar companion and the contamina- 
tor with a template based on the high-resolution spectrum 
of the G2IV star Hyi (H P 2151) from the UVESPOP li- 
brary (Bag nulo et al.ll2003l) . For our fitting, we convolved the 
UVES spectrum with a truncated Gaussian to match the res- 
olution of the observed spectra corresponding to the ~ 0"7 
seeing and 0"7 slit. We ignored rotational and orbital broad- 
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3 From Fig. 1 of |Callanan et aT| fT995l) , we read off R = 19.89 ± 0.03 for 
the contaminator, with a quoted systematic uncertainty of 0.05 mag. This is 
somewhat brighter than what we measured from the HST data. We use the 
fainter magnitude to obtain a conservative limit on the radius. 
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FIG. 2. — Radial- velocity measurements of the companion of PSR 
B 1957+20 as a function of orbital phase. The colored circles mark veloc- 
ity measurements from different nights, with errors as indicated. They are 
repeated in black for clarity. The drawn curve is the best-fit circular orbit, 
and the dotted one the pulsar orbit inferred from radio timing. 

ening, which are well below our resolution for the blue spec- 
tra and just comparable to that of the red spectra (for com- 
panion radius R c ~ 0.25 Rq, one finds vsin/ = 2irR c sini/P ~ 
30kms _1 sin/; for integration time f; nt = 1800s, the maxi- 
mum smearing is Kgb S sin(27rfi nt /f) ~ 80kms _1 , where K \, s = 
324kms _1 is the radial-velocity amplitude [see below]). 
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The template was fitted for a grid of velocities between 
-600 and +600 kms" 1 with a step size of 5 kms" 1 , at each ve- 
locity fitting for the normalization and possible variation with 
wavelength using a quadratic function. Typical reduced \ 2 
values for the best fits were Xred — 1-2. The best fit velocity 
was determined using a quadratic fit to the % 2 values within 
40 kms" 1 of the minimum. Looking at the results for the con- 
taminator, it is clear that systematic variations are present, 
with root-mean-square of 13 kms -1 , likely due to small shifts 
in placement in the slit and/or uncorrected atmospheric dis- 
persion. As our velocities, we thus take the difference be- 
tween the velocities inferred for the companion and the con- 
taminator; these velocities and their corresponding uncertain- 
ties are listed in Table [TJ 

To see how sensitive our results are to our choices, we 
tried fitting only part of the blue spectrum, and using differ- 
ent UVESPOP stars or model atmospheres as templates. We 
found that the velocities were consistent to within ~ la, which 
we will use as an estimate of the associated systematic uncer- 
tainty below. We also tried fitting the red spectra, but found 
that for most, the absence of the expected lines combined with 
the poor signal-to-noise did not allow us to obtain a reliable 
velocity. 

We fitted the velocities with a circular orbit using epoch of 
ascending node T = MJD 48196.0635242 and orbital perio d 
P = 33001.91484s dArzoumanian. Fruchter, & Tavlori[l994T) . 
but found we could obtain a good fit only if we left the 
phase free (Fig. 0. Compared to the prediction, ascending 
node occurs ~ 350s later (phase offset 0.01 1 ± 0.002). How- 
ever, the orbital per i od is k nown to vary quasi-periodically 
dArzoumanian et aT] 119941; iNice. Arzoumanian. & Thorsetj 
2000); indeed, fro m the first and second orbita l period deriva- 
tives measured by lArzoumanian et all d 19941) . the expected 
phase offset is 0.0406. Keeping the phase free, our fit has 
Xred = 0-92 for 12 degrees of freedom (15 measurements, 3 
parameters). The observed radial-velocity amplitude is K Q \, S = 
324 ± 3 kms" 1 and the systemic velocity, measured relative to 
the contaminator, is 7 or)s = 1 10 ±5 kms -1 (here, we multiplied 
the formal errors with v2 to account for the systematic un- 
certainty in the velocities related the choice of template and 
fitting region discussed above). 

The above systemic velocity is relative to the velocity of 
the contaminator. We tried to measure the latter in two dif- 
ferent ways. First, a simple average of the velocity mea- 
surements from the blue spectra yields -39 ±4 kms" 1 (where 
we correcte d for the velocity of 22.7 ± 0.9kms _I of (3 Hyi, 
lEvansI H967L we verified that we obtained the same result 
within <2kms"' using other UVESPOP stars and model at- 
mospheres). Possible evidence that this is not reliable, how- 
ever, comes from the large scatter in the velocities from indi- 
vidual spectra (see above). A similar scatter is seen in the O I 
5577 sky emission line. The latter also shows an average off- 
set, of 6 kms" 1 , but it is not clear whether one can apply this, 
since the line is very close to the red edge of the blue spectra, 
where the wavelength solution may be less reliable. 

As an alternative, therefore, we measured the velocities 
of the contaminator using the red spectra, where, given the 
higher resolution, any shifts due to offsets in the slit and/or 
atmospheric dispersion should be smaller. Here, we unfortu- 
nately could not use the UVESPOP spectra, since these have a 
gap over the Ca II IR triplet, and thus we determined velocities 
relative to a T e ff = 6000 K, logg = 4.5 model spectrum from 



IZwitter et al.l d2004). We find that the scatter in these veloc- 
ities is substantially smaller (root-mean-square of 6kms _1 ), 
and that the (weighted) average velocity is -25.1 ± 1.5 kms" 1 . 
From the OH sky line at 8344.602 A, we infer that any off- 
set due to wavelength calibration errors are small; including 
these, we find an average velocity of -25.6 ± 1.3 kms" 1 . 

As a check on these numbers, we can compare the ve- 
locities with what is expected for a star that is following 
simple Galactic rotation. Assuming a fiat rotation curve 
with = 2 20km s" 1 a nd a distance to the Galactic center 
of 8.5 kpc dCoxl 12000). as well as a peculiar velocity of 
the Sun relative to the lo cal standard of rest of (U ,V,W) = 
(10.00,5.25, 7. 17)kms" 1 (iDehnen & Binneyl [19981) . we find 
that for distances of 6, 8, and lOkpc for the contaminator, the 
expected radial velocities are +16.7, -1, and -24 kms" 1 , con- 
sistent with our measurements. (Arguably, one should con- 
sider asymmetric drift in our estimates, in which case the ve- 
locities would be more negative by ~ 15kms _1 .) 

Overall, we believe the contaminator velocity from the red 
spectra is more reliable. For PSR B1957+20, the implied sys- 
temic velocity is 84 ±5 kms" 1 . 

4. RADIUS AND INCLINATION 

The companion is irradiated by the pulsar and presents a hot 
and a cold side. Hence, the center of light does not coincide 
with the center of mass, but is shifted somewhat towards the 
pulsar. As a result, the spectroscopic observations underes- 
timate the true radial- velocity amplitude. The correction de- 
pends on the stellar radius, the temperature distribution, and 
the inclination. Here, we discuss the available constraints on 
these properties. 

The lightcurve provides strong constraints on the system 
parameters, especially if simplifying assumptions can be 
made for the temperature distribution. For P SR B 1957+20, 
high-q ualit y lightcurves were p resented by ICallanan et al.1 
(1995) and Rey nolds etafl (120071) . These authors also fit their 
lightcurves with models. 

In these lightcurve synthesis models, it is assumed that the 
companion's shape is that of an equipotential surface, and that 
its temperature distribution is given by some background tem- 
perature (modified suitably by gravity darkening) that is in- 
creased by irradiation by an isotropic pulsar wind such that the 
outgoing flux equals the sum of the bac kground and irradia- 
tion flu xes (for a detailed description, see lOrosz & Hauschildtl 
(2000), who wrote the ELC code used by [Rey nolds et al. 
(2007)). It is assumed the irradiation does not affect the tem- 
perature structure of the atmosphere (the "deep heating" ap- 
proximation), such that each surface element can be taken to 
radiate as predicted by a model atmosphere for a single star. 
The main free parameters are the background temperature and 
the irradiating flux (which set the temperature distribution), 
the extent to which the companion fills its Roche lobe (which 
d etermines its shape), an d the inclination of the orbit. 

[Reynolds et afj (120071) find that such models reproduce the 
lightcurves in detail, and they infer that the companion nearly 
fills its Roche lobe, up to a filling factor R m se/RL\ in the range 
0.81 < Rnose/Ru < 0.87 (where R nose is the radius of the star 
in the direction of the pulsar, and the distance to the inner 
Lagrangian point); this corresponds to a volume-equivalent 
radius R 2 in the range 0.946 < Ri/Rrl < 0.974 (where R RL 
is the volume-equivalent radius of the Roche lobe). They also 
infer an inclination i in the range 63° < i < 67°. 

Formally, the above ranges are at the 3a level 
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dRevnolds et al.l 12007b . This, however, does not take 
into account uncertainties in the models or the extent to 
which the underlying assumptions hold. For instance, the 
models may have temperature-dependent missing opacities, 
the metallicity or hydrogen abundance may not be solar, 
the pulsar wind that irradiates the companion may not be 
isotropic (though the lightcurve is remarkably symmetric), or 
the deep heating approximation may not be valid (although 
our blue spectra suggest it is not bad). Given these issues, we 
will treat the ranges as la uncertainties below. 

In addition to the above, an uncertainty that is more diffi- 
cult to constrain relates to the extent to which heat is redis- 
tributed. Any redistribution would lead to a smoother tem- 
perature distribution, and thus it would require a larger in- 
clination to obtain the same observed modulation amplitude. 
In principle, it should be possible to constrain the tempera- 
ture distribution directly by fitting multi-band lightcurves and 
spectr a simultaneously, sim ilar to what has been done for NN 
Ser by Par sons et alj (120101) . which is also strongly irradiated 
(and for which no evidence for heat transport was found). We 
hope to pursue this in the future. 

Here, we will try to use the observations to set limits. For 
our purposes, the most important quantities are the companion 
radius, which determines the correction to the radial-velocity 
amplitude (see below), and the inclination, on which the final 
masses depend as 1 / sin 3 i. 

For the radius, our spectra yield independent clues. First, 
the surface gravity is limited to log g 2 < 4. Since the min- 
imum companion mass is M 2 . m \ n = Q.Q22Mq (from the pul- 
sar mass function and our observed radial-velocity ampli- 
tude), one infers R Xg = (g 2 /GM 2 ) 1/2 > 0.25Rq. Similarly, 
from the distance limit inferred from the reddening, we found 
R 2 ,d>R e n> 0.25 R Q . 

These radii can be compared with the radius of the Roche 
lobe. For a companion much less massive than the pul- 
sar, R RL ~ 0.46a(M 2 /[Mi +M 2 ]) 1/3 dPaczvnskil 119711) . The 
scaling yields the well-known result that the mean density 
of the Roche lobe is determined just by the period; nu- 
merically, p RL = 0.185gcirr 3 (P/ld)- 2 dEggletonlll983l) . For 
PSR B1957+20, 75 RL = 1.27gcm -3 , and hence the size of the 
Roche lobe is fl RL = (3^^/47^) 1/3 = 0.29 R & (where the 
numerical value is for the minimum mass). Thus, we con- 
clude that R 2 . g /RRL > 0.86. This is an overall lower limit, 

since R 2 ^/Rrl oc M^ 6 . Our limit is consistent with what 
was inferred from the lightcurve fit (as well as from the- 
oretica l considerations of the cau se of orbital period varia- 
tions; Applegate & Shaham 1994). It implies a lower limit 
^nose/^n ^ 0.7, which we will use below. 

For the inclination, it is more difficult to set stringent lim- 
its. The fact that at low radio frequencie s, a symmetric eclipse 
is seen that lasts ~ 8% of the orbit (Fruchter et al.l fl990t 
iRyba & Tavlorlll991[). while no eclipse is seen at high fre- 
quencies ((Fruchter & Gossl 19921) . suggests that the conditions 
along the line of sight do not change too strongly during the 
eclipse. This seems easier to understand if the line of sight 
does not pass close to the companion, i.e., if the inclination is 
intermediate, i ~ arccos(/?£/a) ~ 75° (for fractional eclipse 
radius Re/o~ 0.087r), consistent with the inference from the 
lightcurve. Strictly, however, the eclipses only set a weak up- 
per limit: the absence of eclipses at high frequencies implies 
/<arccos(# 2 /a)<85°. 

A lower limit to the inclination can be set from the large 



brightness contrast between su perior and inferior con junction 
(a factor 100 in the R band, IRevnolds et all [2007b . which 
shows that at inferior conjunction at most a small part of the 
irradiated hemisphere is visible. Conservatively, we estimate 
that this requires i > 50°. A similar lower limit seems reason- 
able from the fact that the radio eclipses last long; for lower 
inclinations, it is difficult to envisage a physical eclipse region 
that does not extend all the way to the pulsar. 

Overall, we conclude that a ll observations a re con sistent 
with the model inferences of IRevnolds et alj ([2007), with 
^nose/^Li = 0.84± 0.03 and i = 65 ± 2deg. Considering possi- 
ble systematic uncertainties, a secure constraint on the radius 
seems to be 0.7 < Rnose/Ru < 1, while for the inclination the 
constraint is 50 < i < 85deg. 

5. MOTION OF THE CENTER OF MASS 

The ratio between the observed radial-velocity amplitude 
of the center of light and the actual one of the center of mass 
can be written as K a \, s /K 2 = 1 — / e ff^nose/^2> where the ef- 
fective normalized emission radius f e s is constrained to be 
between (uniform emission) and 1 (emission from the tip 
of the star facing the pulsar only). It should depend pri- 
marily on the surface brightness and line strength distribu- 
tions in the observed band, with minor contributions arising 
from the exact shape of the star (determined by the mass ra- 
tio, filling factor, and the degree of co-rotation) and the or- 
bital inclination. Indeed, for the somewhat similar situation 
of irradiation-indu ced Bowen emission l ines o n a Roche-lobe 
filling companion, Mufio z-Darias et all d2005l) found that the 
"K-correction" factor depends mainly on mass ratio (which 
determines the radius) and is nearly independent of inclina- 
tion. 

We investigated this hypothesis by adapting a lightcurve 
synthesis code to de al with spectra (the code was used by 
Stapp ers et alj d 19991) to model pulsar irrad iation, and is sim- 
ilar to that of Orosz & Hauschildt (2000)). The code pro- 
duces not just fluxes, but also synthetic spectra, by summing 
Doppler-shifted spectra over all surface elements. While we 
do not yet have a suitable set of atmosphere models in hand 
to model our observations reliably, we have used the code to 
estimate the effect on the radial-velocity amplitude, by gen- 
erating mock spectra for different sets of binary parameters. 
We generated spectra at the same orbital phases as our obser- 
vations, determined radial velocities usin g the same method 
as done for the real observations (§ 13. 2b . and fitted circular 
orbits using the same weightsQ 

At a fixed mass ratio, we found that the orbital inclina- 
tion had negligible effect on / e g , < 1% over the range 50 < 
i < 90deg. The effects of the filling factor, for the range 
0.7 < Rnose/Ru < 0.95, were a bit larger, though still small, 
at ~7%. The influence of the precise values of strength of the 
irradiation is ^6%, for the range of irradiation strengths that 
give front-side temperatures consistent with our spectral type 
(6000-6500 K); the backside temperature d oes not matter 
much, since it is constrained to be ~ 2900 K dRevnolds et alJ 
120071) . much too low to contribute B-band flux. As expected, 
with our choice of scaling with R nose , the effect of the mass ra- 
tio is very small, < 1% for the range of values that are able to 
yield the observed K Q \, S . Overall, we infer f e g ~ 0.60 ± 0.04. 

4 For a distorted, irradiated model, the predicted radial-vel ocity curve is 
not n ecessarily circular (see, e.g., the curve for NN Ser of Parsons et all 
120101) . But any resulting systematic effects are corrected for by fitting cir- 
cular orbits to observed and model velocities at the same phases and with the 
same weighting. 
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For comparison, looking from the side at a spherical star 
that is dark on one side and emits isotropically at a uni- 
form temperature on the other, simple integration yields / e ff = 
4/3-7T = 0.42. Assuming instead a temperature distribution 
T = 7o cos 1 / 4 9, as expected for irradiation by a parallel beam, 
we find / e ff = 0.62 for 7q = 6400 K (assuming black-body 
emission and a linear limb darkening law with u = 0.6). This 
is very similar to what we find from our model, showing that 
the increase in brightness dominates over effects such as an 
increase in line strength with decreasing temperature. 

Before using the above to estimate K2, it is useful first 
to consider the limits. For uniform emission, K2 = = 
324 ± 3kms _I , which yields a lower limit to the mass ratio 
for the system q mm = M\ /M2 = K \, s /K\ = 63.6 ± 0.6 (where 
the pulsar's radial- velocity amplitude K\ = 2na\ sin//P or b = 
5.09272 ± 0.00004kms~ 1 is known from radio timing; 
lArzoumanian et aD fl994). At the other extreme, for light 
emitted by the nose of a Roche-lobe filling star, i.e., / e ffR n ose = 
R LU one finds K obs = 0.843# 2 , i.e., K 2 = 384±4kms _1 and 
<7max = 75.4 ±0.7. We will take these values as limits in §|6] 
(accounting for the fact that neither can be close to realistic 
by ignoring the uncertainty on K b s ). 

For our more general case, combining our estimate of / e ff = 
0.60 ±0.04 above with R nose /R L i = 0.84 ±0.03, and solv- 
ing for K2 (taking into account that Ru/a depends on q), 
one finds K2 = 353 ±4kms _1 . For the conservative range in 
companion radius, 0.7 < ^ n0S e/^Li < 1, we find 348 ±4 < 
K2 < 358 ±4kms _I , i.e., it corresponds to a 5kms _1 uncer- 
tainty in K2- Adding this in quadrature to the 2a uncer- 
tainty of 8kms _I arising from the uncertainties in K Q \, S and 
/ e ft, we infer a conservative range in radial-velocity ampli- 
tude 343 < K2 < 363 km s" 1 . From the above, we conclude 
that the mass ratio is q = 69.2 ± 0.8, and that a conservative 
range is 67.3 < q < 71.3. 

6. INFERRED MASSES AND CONCLUSIONS 

In Fig. [3] we show our constraints on the masses. One 
sees that most likely, PSR B 1957+20 is mass ive, with Mpsr = 
2.40M Q . Taking the inferences from iReynolds et al.l (2007) 
on the lightcurve and the corresponding inclination, and our 
correction to the radial-velocity amplitude at face value, the 
formal uncertainty is small, ^0.12M Q . 

As discussed in §|4] however, the lightcurve modelling re- 
lies on a number of assumptions, especially that there is no 
heat transport over the face of the star. From our conservative 
constraints on both the inclination and the mass ratio, we find 
a lower limit to the mass of 1 .66M Q . 

Thus, from our work we conclude that PSR B 1957+20 cer- 
tainly is more massive than the canonical 1.35M Q , and likely 
substantially more massive. Indeed, it may w ell be more mas- 
sive e ven than PSR J1614-2230 ( 1 .97 ± 0.04: lDemorest eFall 
2010), and thus allow even more stringent constraints on the 
equation of state. The large mass also suggests a large amount 
of mass was transferred in the preceding phase as an X-ray 
binary, although this conclusion depends on the initial mass. 
However, even if tha t were as high as the mass found for 
Vela X-l (~ 1.9M Q , iBarziv et al.ll200Th . our measurements 
suggest the pulsar has accreted about half a solar mass0 

To confirm the high inferred mass will require more secure 
constraints on the orbital inclination and, to a lesser extent, 

5 Since Vela X-l has a massive companion, its mass must still be close to 
the one it formed with. 
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FIG. 3. — Mass-mass diagram for PSR B1957+20 and its companion. The 
cross indicates our best-fit solution with la uncertainties, and the surround- 
ing blue parallelogram the convervative region including our best estimate 
of possible systematic uncertainties (see text). The physically allowed re- 
gion (light blue) is limited by the constraints sin (' < 1 (thick, blue line), and 
< i^obs < K2 (thick, red lines). Here, the second constraint arises because 
the observed radial-velocity amplitude K2 is measured using light emitted 
from the side facing the pulsar, and the center of light cannot be further away 
from the pulsar than the center of mass, or closer to the pulsar than the first 
Lagrangian point LI (with velocity amplitudes K2 and ATli, respectively). 
For reference, we show contours of constant inclination i (dotted), calculated 
using the pulsar mass function. 

the correction factor for the radial-velocity amplitude. For this 
purpose, most important would be to model the lightcurve in 
more bands, and to check explicitly what inclinations are pos- 
sible for less-constrained, perhaps even arbitrary temperature 
distributions. It may be especially valuable to model the spec- 
tra at the same time, thus avoiding the indirect calculation of 
a correction factor. We are currently working on rewriting our 
code for this purpose. 

An improved, nearly model-independent constraint on the 
mass ratio could be obtained from a near-infrared radial- 
velocity curve. Since in the near-infrared the contribution of 
the cold side to the light budget is more important, the radial- 
velocity amplitude of the center of light at infrared wave- 
lengths should be much closer to that of the companion's cen- 
ter of mass, and hence the uncertainty in the corrections much 
less important. 

Furthermore, one could determine the projected rotational 
velocity vsim from high-resolution spectra. This would allow 
one to check the predictions from the models, in particular for 
the mass ratio and the filling factor, on which vsim depends 
most strongly. 

Finally, radio observations could provide a complemen- 
tary improved constraint on the inclination, from mapping 
the eclipse at a larger ran ge of frequencies than was done 
by lFruchter & Gossl (119921) . and making use of the large in- 
creases in sensitivity, especially at high frequency, that have 
been made over the last decades. 

The data presented herein were obtained at the W.M. Keck 
Observatory, which is operated as a scientific partnership 
among the California Institute of Technology, the University 
of California and the National Aeronautics and Space Admin- 
istration. The Observatory was made possible by the generous 
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financial support of the W.M. Keck Foundation. We also used 
data from the UVES Paranal Observatory Project UVESPOP 
(ESO DDT Program ID 266.D-5655). We made extensive use 



of SIMBAD and ADS. 

Facilities: Keck (LRIS) 
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